# distance to metro stations, parks, schools, and hospitals in bcn
# -------------------------------------------------------------------- #
# prelims ####
# -------------------------------------------------------------------- #
# clear environment
rm(list = ls())

# import relevant libraries
library(haven)
library(rgdal)
library(raster)
library(readxl) 
library(classInt)
library(tmap) 
library(Hmisc) 
library(GISTools) 
library(tidyverse)
library(units)
library(sf) 
library(RANN) 
library(fields) 
library(geosphere) 

# paths
# main directory
dir <- paste0("PATH_TO_MAIN_DIRECTORY_HERE/")
setwd(dir)
# ajuntament de barcelona
SecCenPath <- paste0(dir,"orig/aj_barcelona/0301040100_BCN_UNITATS_ADM/0301040100_SecCens_ADM_ETRS89.shp")
# generalitat
GenPath <- paste0(dir,"orig/generalitat/XT_FERRO_ESTACIONS/XT_FERRO_ESTACIONS.shp")
# vars to store
varscensus <- c('DISTRICTE','BARRI','AEB','SEC_CENS','GRANBARRI','ZUA','LITERAL','AREA')

# -------------------------------------------------------------------- #
# metro stations ####
# -------------------------------------------------------------------- #

# barcelona shape
GisDataSecCen <- st_read(dsn = SecCenPath) %>%
  select(c(all_of(varscensus)))
CRS <- st_crs(GisDataSecCen)

# fgc stataions
CatStations <- st_read(dsn = GenPath) %>%
  st_transform(CRS)  

# tmb stations
TMB <- read_dta(paste0(dir, "data/int/metro_stations.dta")) %>%
  st_as_sf(coords = c("lon", "lat"), crs = '+proj=longlat +datum=WGS84') %>%
  st_transform(CRS)  

# keep stations in barcelona only
BCNStations <- st_intersection(CatStations, GisDataSecCen)

# get distance to closest stations from centroid of each census tract
Cent <- st_centroid(GisDataSecCen) %>%
  select(c("geometry"))

# calculate distances
NTracts <- nrow(Cent)
NStat <- nrow(BCNStations)

# matrix of distances
Distances <- data.frame()

# fill matrix
for(stat in 1:NStat){
  
  for(tract in 1:NTracts){ 
    # calculate distance between tract and station
    Distances[tract,stat] <- st_distance(Cent[tract,"geometry"],BCNStations[stat,"geometry"])
  } 
  
  # rename column
  colname <- paste("V",stat, sep="")
  statname <- paste(BCNStations$CODIPOL[stat], sep="") %>%
    str_replace_all("[[:punct:]]", "")
  names(Distances)[names(Distances)==colname] <- statname
  
} 

# put the two datasets together
GisDataSecCen$ID = seq.int(nrow(GisDataSecCen)) 
# remove geometry
GisDataSecCen$geometry <- NULL

# merge
Distances$ID = seq.int(nrow(Distances)) 
CTractDist <- merge(GisDataSecCen, Distances, by="ID")

# correct variable names
names(CTractDist) <- gsub("-|\\.|\\/|'|\\[|\\]", "", names(CTractDist))
names(CTractDist) <- gsub(" ", "", names(CTractDist))

# store data
filename_R <- paste0(dir,"data/int/distance_stations_tract.RData")
filename_dta <- paste0(dir,"data/int/distance_stations_tract.dta")
write_dta(data = CTractDist, filename_dta)
save(CTractDist, file = filename_R)

# remove
rm(CatStations, TMB, BCNStations, CTractDist, Cent, Distances, GisDataSecCen)

# -------------------------------------------------------------------- #
# schools ####
# -------------------------------------------------------------------- #

# barcelona shape
GisDataSecCen <- st_read(dsn = SecCenPath) %>%
  select(c(all_of(varscensus)))
CRS <- st_crs(GisDataSecCen)

# schools
Schools <- read_dta(paste0(dir, "data/int/bcn_schools.dta")) %>%
  # keep only actual schools
  subset(school == 1)  %>%
  # transform to spatial
  st_as_sf(coords = c("lon", "lat"), crs = '+proj=longlat +datum=WGS84') %>%
  # transform coordinateS
  st_transform(CRS)  

# keep schools in barcelona only
BCNSchools <- st_intersection(Schools, GisDataSecCen)

# get distance
Cent <- st_centroid(GisDataSecCen) %>%
  select(c("geometry"))

# calculate distances
NTracts <- nrow(Cent)
NStat <- nrow(BCNSchools)

# matrix of distances
Distances <- data.frame(matrix(nrow=NTracts, ncol=NStat))

# fill matrix
for(stat in 1:NStat){

  for(tract in 1:NTracts){ 
    # calculate distance between tract and station
    Distances[tract,stat] <- st_distance(Cent[tract,"geometry"],BCNSchools[stat,"geometry"])
  } 
  
  # rename column
  schlID <- paste0(BCNSchools$register_id[stat]) %>%
    # remove special characters
    str_replace_all("[[:punct:]]", "")
  # rename
  names(Distances)[stat] <- paste0("SCL",schlID)
  
} 

# put the two datasets together
GisDataSecCen$ID = seq.int(nrow(GisDataSecCen)) 
# remove geometry
GisDataSecCen$geometry <- NULL

# merge
Distances$ID = seq.int(nrow(Distances)) 
CTractDist <- merge(GisDataSecCen, Distances, by="ID")

# correct variable names
names(CTractDist) <- gsub("-|\\.|\\/|'|\\[|\\]", "", names(CTractDist))
names(CTractDist) <- gsub(" ", "", names(CTractDist))

# store data
filename_R <- paste0(dir,"data/int/distance_schools_tract.RData")
filename_dta <- paste0(dir,"data/int/distance_schools_tract.dta")
write_dta(data = CTractDist, filename_dta)
save(CTractDist, file = filename_R)

# remove objects
rm(Schools, BCNSchools, CTractDist, Cent, Distances, GisDataSecCen)

# -------------------------------------------------------------------- #
# parks ####
# -------------------------------------------------------------------- #

# barcelona shape
GisDataSecCen <- st_read(dsn = SecCenPath) %>%
  select(c(all_of(varscensus)))
CRS <- st_crs(GisDataSecCen)

# parks
Parks <- read_dta(paste0(dir, "data/int/bcn_parks.dta")) %>%
  st_as_sf(coords = c("lon", "lat"), crs = '+proj=longlat +datum=WGS84') %>%
  st_transform(CRS)  

# keep parks in barcelona only
BCNParks <- st_intersection(Parks, GisDataSecCen)

# get distance to closest park from centroid of each census tract
Cent <- st_centroid(GisDataSecCen) %>%
  select(c("geometry"))

# calculate distances
NTracts <- nrow(Cent)
NStat <- nrow(BCNParks)

# matrix of distances
Distances <- data.frame(matrix(nrow=NTracts, ncol=NStat))

# fill matrix
for(stat in 1:NStat){

  for(tract in 1:NTracts){ 
    Distances[tract,stat] <- st_distance(Cent[tract,"geometry"],BCNParks[stat,"geometry"])
  } 
  
  # rename column
  parkID <- paste0(BCNParks$register_id[stat]) %>%
    # remove special characters
    str_replace_all("[[:punct:]]", "")
  # rename
  names(Distances)[stat] <- paste0("PRK",parkID)
  
} 

# put the two datasets together
GisDataSecCen$ID = seq.int(nrow(GisDataSecCen)) 
# remove geometry
GisDataSecCen$geometry <- NULL

# merge
Distances$ID = seq.int(nrow(Distances)) 
CTractDist <- merge(GisDataSecCen, Distances, by="ID")

# correct variable names
names(CTractDist) <- gsub("-|\\.|\\/|'|\\[|\\]", "", names(CTractDist))
names(CTractDist) <- gsub(" ", "", names(CTractDist))

# store data
filename_R <- paste0(dir,"data/int/distance_parks_tract.RData")
filename_dta <- paste0(dir,"data/int/distance_parks_tract.dta")
write_dta(data = CTractDist, filename_dta)
save(CTractDist, file = filename_R)

# hospitals ####
# -------------------------------------------------------------------- #

# barcelona shape
GisDataSecCen <- st_read(dsn = SecCenPath) %>%
  select(c(all_of(varscensus)))
CRS <- st_crs(GisDataSecCen)

# hospitals
Hospitals <- read_dta(paste0(dir, "data/int/bcn_hospitals.dta")) %>%
  st_as_sf(coords = c("lon", "lat"), crs = '+proj=longlat +datum=WGS84') %>%
  st_transform(CRS)  

# keep hospitals in barcelona only
BCNHosp <- st_intersection(Hospitals, GisDataSecCen)

# get distance to closest hospital from centroid of each census tract
Cent <- st_centroid(GisDataSecCen) %>%
  select(c("geometry"))

# calculate distances
NTracts <- nrow(Cent)
NStat <- nrow(BCNHosp)

# matrix of distances
Distances <- data.frame(matrix(nrow=NTracts, ncol=NStat))

# fill matrix
for(stat in 1:NStat){

  for(tract in 1:NTracts){ 
    # calculate distance between tract and station
    Distances[tract,stat] <- st_distance(Cent[tract,"geometry"],BCNHosp[stat,"geometry"])
  } 
  
  # rename column
  hospID <- paste0(BCNHosp$register_id[stat]) %>%
    # remove special characters
    str_replace_all("[[:punct:]]", "")
  # rename
  names(Distances)[stat] <- paste0("HSP",hospID)
  
} 

# put the two datasets together
GisDataSecCen$ID = seq.int(nrow(GisDataSecCen)) 
# remove geometry
GisDataSecCen$geometry <- NULL

# merge
Distances$ID = seq.int(nrow(Distances)) 
CTractDist <- merge(GisDataSecCen, Distances, by="ID")

# correct variable names
names(CTractDist) <- gsub("-|\\.|\\/|'|\\[|\\]", "", names(CTractDist))
names(CTractDist) <- gsub(" ", "", names(CTractDist))

# store data
filename_R <- paste0(dir,"data/int/distance_hospitals_tract.RData")
filename_dta <- paste0(dir,"data/int/distance_hospitals_tract.dta")
write_dta(data = CTractDist, filename_dta)
save(CTractDist, file = filename_R)

# -------------------------------------------------------------------- #
# closing ####
# -------------------------------------------------------------------- #

rm(list = ls())
